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The development of the SHELPC finite element computer program is 
detailed. This program is specialized to simulate the nonlinear material 
behavior which results from combustor liner -hot streaks". This problem 
produces a nonlinear Fourier Series type loading on an axisynmetric 
structure. Example cases are presented. 


One problem which is unique to Aircraft Gas Turbine Engines (AGTE’s) is 
the thermal hot streak problem experienced by certain combustor liners. These 

pressure shells whose function is to contain and 
promote the combustion process. Because of the high temperatures involved in 

functioning of these liners depends on their 
‘>yJo'^er temperature air from the compressor. To complicate 
nnt S® temperatures generated by the combustion process are 

liner?^^ circumferential directions of the combustor 

The temperature variation in the axial direction is the variation normally 
shown in textbooks and calculated from the thermodynamic laws based on 
pressure/yolume considerations. The circumferential variation is a more 
SS dimensional mixing problem. It is a function of the number of 

fuel nozzles employed and the geometry and thermodynamics of the combustor. 

^e result is a circumferential variation in the metal temperature of these 
liners which can be approximated with Fourier Series. 

lives of these liners are primarily dictated by the 
l^sponse to these temperature variations, which produce 
ther^l stresses. The stresses produced by the differential pressures across 

the nrlhiJH*^^! to small values and are secondary contributors to 

Because of the high temperatures and large temperature 

J^® response problem is nonlinear and time 

^®®” active in developing internal computer tools for 
twft types of unique AGTE problems. Our capabilities include both 

two-dimensional and three-dimensional nonlinear finite element computer 
programs. 
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The thermal stress problem generated by the axial variation In metal 
temperatures can be handled quite readily by a two*d1mens1onal (axlsyranetrlc) 
analysis. However, the thermal stresses generated by the circumferential 
variation In metal temperatures Is truly a three-dimensional problem. As 
such, one method of attack would be to develop a simulation using 
three-dimensional finite elements, such as our 8, 16 or 20 noded Isoparametric 
elements. 

This method of attack has drawbacks In economy, accuracy, and utility. 

The number of hot streaks can be quite large (In the 60 's) and the hot to cold 
temperature variations can be In the hundreds of degrees Fahrenheit. This 
would require a large number of 3D elements to produce an acceptable 
simulation. This could also produce bad element aspect ratios, resulting In 
solution accuracy problems. Varying these models In an Iterative design 
process would be manpower and computer Intensive. For these reasons we sought 
other methods of attacking this problem. 

Since the circumferential temperature variations can be represented as 
Fourier Series expansions, this suggests attacking the structural problem by 
means of Fourier Series. This method of attack for linear elastic problems in 
a finite element format was presented by E.L. Wilson In Reference ( j). For 
linear elastic problems the method of superposition of Fourier Series loading 
Is mathematically exact and precise. However the problem In question Is 
nonlinear. 

In Reference (2), E.A. Witmer and J.J. Kotanchik presented a nonlinear 
solution procedure for axisymmetric shells under asyrnnetrlcal loading. In 
their approach, they utilized an Initial strain method with Besseling's 
Isothermal constitutive model In a sublayer format. This then presented a 
method of attacking nonlinear Fourier Series loading problems but with one 
tremendous deficiency; our problem under consideration Is totally 
non Isothermal In nature. 

Subsequently, one of the present authors, R.L. McKnIght, developed a 
variable temperature version of Besseling's constitutive model In a subvolume 
format as presented In References (3) and (4). This has been the primary 
basis of AEBG's internal nonlinear computer tools since 1975. To this 
time-independent plasticity formulation, we added classical creep capability 
with time-hardening, strain-hardening, and life-fraction rules for attacking 
quasi-static time-dependent nonlinear problems. Over the years, through large 
amounts of production usage, these tools have been provided with much 
verification, validation, and operational experience for both 2D and 3D 
problems. With this background, we decided to reattack the combustor liner 
hot streak problem: 


SHELPC 

The following Is a synopsis of the theoretical development of the 
combustor liner program, SHELPC. A complete development will be found In the 
Ph.D. Dissertation of P.C. Chen, to be published later. 

As in Reference (1), the finite element employed Is the triangular ring 
element shown in Figure 1. This is an extension of that element commonly used 
for plane stress, plane strain, and axisymmetric analysis. In those cases, 
there are two degrees of freedom, or displacements at each node. To allow for 
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the circumferential variation in load a third degree of freedom, the 
circumferential or hoop displacement, is introduced. These three displacement 
components are assumed to be linear functions of position in the R-Z plane. 
This is expressed mathematically as follows 


Ur - bi + b2R + b3Z 
Uz = b4 + bgR + bgZ 
Ug ■ by + bgR + bgZ 


( 1 ) 


The strait) components in terms of these displacements are (for small displace- 
ment theory). 
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The loading is now allowed to varying in the circumferential direction in a 
sine or cosine fashion as defined in Table 1. 


Equations (1), then become 


“r ’ 

“z ’ C’4n * '’Sn'* * •’6n^> (3) 

“e = ("Zn * ‘>8n'* * 

“r ' <»ln * ‘>2n“ * »'3n« ^1" Ne 

'*2 " ^>>40 * '•sn'' * ‘>6n^> «) 

* C’zn * ■’sn'* * 
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TABLE 1 


FOURIER SERIES LOADING AND DISPLACEMENT DESCRIPTION 


Syirmetric about a Ant1-synmetric about a 

plane containing the plane containing the 

axis of revolution axis of revolution 


T = 2T^(r,z) cos ne 

Sf = lS^j^(r,z) cos ne 

' iS^^Cr.z) cos ne 

Sq » rSgj^(r,z) sin ne 

= iU^(r,z) cos ne 
r m 

= rU2^( r,z) cos ne 
Ue “ £Ugj^(r,z) sin ne 


T = zT^(r,z) sin n e 
Sp = zS^^(r.z) sin ne 
* zS^j^Cr.z) sin ne 
S9 = ZS9j^(r,z) cos ne 
Ur = zUr^(r,z) sin ne 

Uz = 

Ut * zU^^(r,z) cos ne 


n is the harmonic number. 
T is the temperature. 

S's are the loads. 

U's are the displacements. 
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4 . 1 . Jo nonlinear material behavior, we make the classical assumption 

that the total strains consist of a summation of elastic, plastic, creep, and 
thermal components 

e P c T 
e ■ c + e + e + e 

We also assume that stress Is linearly related to the elastic strain only 

- • r -6 


( 5 ) 


o » C(e - e** - - e^) 

Where C is the elastic Hook's Law matrix. 

The sum of the plastic, creep, and thermal strains are considered as 
Initial strains. 


( 6 ) 

(7) 


. cP + cC + cT 


or 


C(e - e^) 


( 8 ) 

(9) 


Applying the precepts of Reference (2). we make use of the Principle of 
Stationary Total Potential Energy. ^ 

The strain energy density Is expressed as 

0 - % J c* 


or 


D » Js[C(c - c^)^3 (c - ch 

0 • >5 (e - C(e - c^) 

The total potential energy of one ring element Is 
• /y Me - C (c - e^)dV - 

^ere 6 Is the generalized displacements and F Is the generalized external 
forces. 


( 10 ) 

( 11 ) 

( 12 ) 

(13) 


For a system of ring elements, this development leads to 
Kfi « F + f1 

where 


( 14 ) 


K ■ elastic stiffness matrix for 
discretized system 
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6 « generalized displacements of the 

discretized system 

F ■ applied generalized forces 

• generalized forces due to Initial 
strains 

Now, to Introduce the Fourier Series loading Into this system of 
equations. Our assumption Is that the total strain can be considered to be 
made up of the sum of a certain number of A*ser1es and B~ser1es Fourier 
components. The A*ser1es components are those symnetrlcal about a plane 
containing the axis of rotation. The B*ser1es components are those 
anti symmetrical about a plane containing the axis of rotation (see Table 1). 
Thus we assume 


U(R,2.0) - r uj5 (R.Z) cos Ne 
+ I uj (R.Z) Sin Ne 

and therefore 

c(R.Z,e) • z cjj (R.Z) cos Ne 
+Z ej (R.Z) Sin Ne 

or 

c(R.Z.e) « c\R,Z,e) + c®(R,Z.e) 


and, using Equation 

(5) 


c(R»Z,e) * 

c®^R.Z,0) 

+ e^^(R,Z,e) + c^^(R,Z,e) 

+ 

e^^R.Z.e) 

+ c®®(R,Z,e) 

+ 

c^®(R,Z.e) 

+ e“(R,Z.e) 

+ 

c*^(R,Z,0) 



(15) 


(16) 


(17) 


(18) 


For a discretized model, the system of equations are set up and solved a 
harmonic and a series at a time 



B 

N 




IB 


( 20 ) 
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( 21 ) 


The nodal displacements are then given by 

6(R,Z.e) ■ I cos Ne + z 6® Sin Ne 


NONLINEAR SOLUTION SCHEME 


A given axlsynmetrlc geometry with thermal and mechanical loadings which 
can be simulated by Fourier Series Is first discretized by triangular rings. 
The pertinent system of equations, (19) and (20), are set up and solved, first 
assuming elastic behavior. This gives a first approximation for the nodal 
point displacements. With these results, which are the amplitudes of Fourier 
Series displacement components^ we can now determine the displacements at any 
circumferential location. 


A minimum period Is selected based on the lowest order harmonic other than 
0 (for the combustor liner this Is determined by the number of fuel nozzles). 
This minimum period Is then approximated by 10 points In the 0-direction. At 
each of these 10 points, each rings nodal point displacements can be 
determined. From this the total strain can be determined for each location. 
Then using the constitutive models as covered In References (3) and (4), an 
initial estimate is made for the Inelastic strains and from these Inelastic 
pseudo-forces, pP and F^. These Inelastic pseudo-forces are then 
approximated by the Fourier Series harmonics. 


F^(R,Z,e) » (R,Z) Cos Ne + (R.Z) Sin Ne 


N 


F^(R,Z,e) * EF,^^^ (R,Z) Cos Ne + EF,^'^ (R,Z) Sin Ne 

The amplitudes of these harmonic pseudo-forces are then added to the Initial 

force vector. ^ ^ TA/« <. .. .. TB/« .. 

Ffj (R,Z,e) « Fj^'”(R,Z) Cos Ne + Fj^ (R,Z) Sin Ne 


CB 


( 22 ) 


+ Fj^**^R,Z) Cos Ne + Fjj*’®(R.Z) Sin Ne (23) 

+ Ff^^(R.Z) Cos Ne + Ff^®(R.Z) Sin Ne 


This new value for the initial force vector Is used to obtain a new series of 
solutions. From these, new Inelastic strains are predicted and from these a 
third Initial force vector. This process Is continued until convergence 
occurs. 


EXAMPLE CASES 

In order to verify the correctness of the solution scheme as outlined, 
several comparisons with finite element methods were made. Two problems 
analyzed are presented, one representing combined thermoplasticity and creep 
under cyclic loading, and the other a verification of the elastic Fourier 
series analysis capability. These cases represent a check on the two basic 
aspects of the combined elastic/plastic technique for harTOnic loading j 
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In the first case, which verifies the cyclic plasticity and creep 
capability, a thick-walled cylinder is subjected to a time-varying pressure 
and temperature loading history (Figure 2). These loadings are the 
predominant types experienced by combustors. Plane strain conditions were 
assumed, the finite element model being composed of triangular ring elements 
as shown in Figure 3. The model was also run using the CYANIDE 2-D computer 
code, subjected to the same temperature and pressure history as shown in 
Figure 2. The comparison between SIELPC and CYANIDE 2-D under these 
axisynnietric loadings was then made, the results being presented in Figures 2 
and 4. Figure 2 shows a comparison of radial displacement versus time for the 
two methods. Figure 5 shows the residual stress distribution through the 
wall at t « 10 hours when the pressure was equal to zero. As can be seen from 
these results, correlation was very good. 

A second test case was used to check the Fourier series method under a 
thermal loading represented by the combination of temperature harmonics of the 
form: 


T(0) - 70 + 30 cos (48) 

In an actual component, such as a combustor, the harmonics used would be a 
function of the number of "hot streaks" around the circumference. The model 
used was again a thick-walled cylinder as shown in Figure 3. To verify the 
response of the structure to this type of loading, the same model was run 
using the CLASS/MASS code which also has harmonic loading capability but 
limited to linear elastic behavior. This problem was run elastically with 
SHELPC in order to achieve a direct comparison. Results are presented in 
tabular form to demonstrate the closeness of the numerical comparison. In 
Table 2, the displacement components from SHELPC and CLASS/MASS are compared 
for the harmonic loading: 

T(9) - 30 cos (49) 

Table 2: Nodal Displacements at 9 ■ 0^ Due to a 
Single Harmonic Thermal Load 
T(0) « 30 cos (40). 


Node 

SHELPC (10”5 inch) 

CLASS/MASS (10"5 inch) 

No. 

*R 


*• 


^2 

6e 

D 

jgm 

-7.726 

8.633 


-7.736 

8.641 

mm 

lli 

7.127 

8.633 

mm 


8.641 

H 

-0.2380 

0 

IQI 

-0.2A45 

Bi 

3.646 

Hi 

7.718 

-8.767 

IBS 

7.719 


9.678 

B 

7.718 

8.767 

9.678 

7.719 


9.678 


In Table 3, the individual harmonic results from SHELPC for the displacements, 
stresses and strains are shown. The results from a combined harmonic run are 
also tabulated. As can be seen from this table, the results are correctly 
superimposed. They also correlate well with the corresponding CLASS/MASS 
values. 


318 












0.420 

0.420 

0.420 

0 

0.174 

0.t)S 

0.187 

-0.0480 

0.594 

0 $ 7 } 

0 607 

0.0480 












0.420 

0.420 

0.420 

0 

0.167 

0.16) 

0.192 

-0.0253 

0.557 

0.585 

0.612 

0.0253 

0.420 

0.420 

0.420 

0 

0.167 

0.16) 

0.192 

1 

-0.025A 

0.587 

0.585 

0.612 

0.0254 

0.420 

0.420 

0.420 

0 

0.1)7 

0.17) 

O.I9> 1 

-0.0320 

0.579 

0.595 

0.615 

-0.0320 


Clc««nt Strcaaet (lO^ p«0 
r 


Et«. 

No. 

■■ 

■■ 

eoh 

■Hi 

1 

Combined 1 

El 

IE] 

IBl 

sa 

«r 

«2 


®lr 

«r 

«2 

«• 

«lr 

1 

D 

0 

0 

0 

-0.549 

-1.00 

-0.244 

0.553 

-0.549 

-1.00 

-0.244 

0.5)3 

2 

0 

D 

0 

0 

-2.592 

-0.631 

-0.005 

-0.292 

-0.592 

-0.631 

-0.005 

-0.92 

m 

0 

0 

0 

B 

-0.591 

-0.631 

-0.005 

-0.293 

-0.591 

-0.631 

-0.005 

-0.93 

B 

0 

0 

a 

B 

-0.664 

-0.289 

0.165 

-0.370 

-0.664 

-0.289 

-0.165 

-0.370 
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Then three more complex nonlinear comparison cases were run to demonstrate 
the capabilities of SHELPC. These cases demonstrate elastic/plastic 
comparisons for axlsyranetrlcally loaded structures as well as harm^lc 
sumnatlon comparisons In the elastic regime. Since no other code which 
the harmonic analysis was available, a direct comparison of nonlinear harmonic 
analysis could not be made. This comparison will rely primarily on 
experimental results from actual combustors. 

To compare the elastic/plastic analysis capability for axlsynmetric 
loadlno. we modeled a typical axlsynmetric combustor lip, as shown In 
Floure 5. It was modeled using 168 3-node finite elements In both SHELPC and 
CYANIDE 2-D programs. The axlsynmetric loading conditions 
compared satisfactorily. Figures 6 and 7 show comparisons of the effective 
stress and effective piastre strain, respectively, as predicted by the two 
orograms. The plot Is for the Inner surface of the cooling louver around the 
edge from element No. 56 to element No. 136 as shown In Figure 8. As can be 
seen from these results, the comparison between the two programs is 
excellent. This case was one of those used to validate the analytical 
development of SHELPC where a direct comparison with axlsynmetric loads Uero 
harmonic) could be made. 

Our second example demonstrates the Fourier series capability of SHELPC 
under a mechanical loading situation. Figure 9 gives the particulars of the 
oroblem. a ring under a harmonic pressure distribution applied at the outer 
diameter. This problem was solved both elastically and plastically for the 
stress-strain properties shown. Figure 10 shows the SHELPC predic^ons for 
the radial displacement at the outer diameter for the ^o cases. The elastic 
results were correlated with hand calculations and CLASS/MASS results. 

As a final validation example, we extended the previous comparison with 
CLASS/MASS to the nonlinear regime. 

Figure 11 shows the problem conditions, an increased thermal loading, and 
the assumed stress-strain curve. The elastic solution was correlated with 
CLASS/MASS. The plastic solution shows the difference In effect between the 
pressure load of the previous example and a thermal hot streak. Figure 12. in 
this case, plastic flow reduces the radial displacement. 

This program is continuing under development with more verification and 
validation cases being pursued. 
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Figure 1. Finite Element Node Ordering. 
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E • 30 X 10 psl 
V ■ 0.3 

a - 6.0 X lO"^ in/(in - 


Node 

No. 

1 

2 

3 

4 

5 


1.0 0.0 

1.0 1.0 

1.5 0.5 

2.0 0.0 

2.0 1.0 
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Effective Stress 


Figure 5. SHELPC Model for Combustor Lip. 



Z Coordinate 

Figure 6. Effective Stress Comparison for Combustor Model. 
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